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Abstract 



r^ ' In order to process a potential moment sequence by the entropy optimization 

"t^ . method one has to be assured that the original measure is absolutely contin- 

uous with respect to Lebesgue measure. We propose a non-linear exponential 
transform of the moment sequence of any measure, including singular ones, so 
that the entropy optimization method can still be used in the reconstruction 
Cn , or approximation of the original. The Cauchy transform in one variable, used 

J^ ' for this very purpose in a classical context by A. A. Markov and followers, is 

-y-v . replaced in higher dimensions by the Fantappie transform. Several algorithms 

p^ ' for reconstruction from moments are sketched, while we intend to provide the 

CO ■ numerical experiments and computational aspects in a subsequent article. The 

^■f-N ' essentials of complex analysis, harmonic analysis, and entropy optimization are 

(^ . recalled in some detail, with the goal of making the main results more accessible 

O^ ' to non-expert readers. 
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1. Introduction 

In sciences and engineering, a particular inverse problem arises often, requir- 
ing approximation of a measure by a density function from knowledge of linear 
data, e.g., integrals of a function basis against a measure. Classical moment 
problem considers integrals of monomials, the power moments, as the set of 
known measurements, while the generalized moment problems extend the ad- 
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missible inputs to integrals of orthogonal polynomials, Fourier basis, wavelets, 
or other functional bases. 

The list of applications of the moment problem is long, ranging from engi- 
neering, through physics, statistics, well into applied mathematics. While pure 
mathematical settings allow for infinite moment sequences, leading to classical 
moment problems of HausdorfF, Hamburger, and Stieltjes, the applied settings 
almost exclusively assume knowledge of only a finite number of moments, which 
is known as the truncated moment problem. 

In early 1980s, statistical physics and signal processing communities rec- 
ognized that a practical solution to the truncated moment problem, which is 
mathematically under-determined, can be found through optimization of the 
Shannon entropy, a nonlinear functional acting on the density of the measure 
[1, 2]. Initial success, in the numerically unfavorable setting of power moments, 
generated sufficient interest to improve on the original method [3-6] and ar- 
rive at a routinely-used method not only in physics, but also in statistics and 
control theory [7-9] . Furthermore, optimization of entropy has been shown to 
be of theoretical importance: it can be used to fully characterize the moment 
sequences representable by densities based on truncated moment data [10] , and 
arbitrarily incomplete moment data [11]. 

However, not every moment sequence is a suitable input for the entropy 
optimization. In particular, it is easy to demonstrate that the moment sequence 
of the Dirac-(S distribution is not a feasible input, as the optimization does not 
converge in that case. Such singular measures captured our focus, as we were 
motivated by potential applications to inverse problems in dynamical systems. 

Measures invariant under evolution of dynamical systems are of particular 
interest, with increasing activity driven by applied problems. On chaotic attrac- 
tors, trajectories of dynamical systems are known to be non-robust to any errors 
and behavior is more reliably represented using statistical methods [12]. Sur- 
prisingly, even in chaotic regimes, the moment data of invariant measures can 
be reliably computed from simulated and experimental trajectories by averaging 
moment functions along them, despite the errors inherent to those procedures 
[13]. Singular invariant measures abound in dynamical systems, e.g., a system 
with an attracting fixed point preserves a Dirac-(5 distribution, whose moments 
are easily computed by averaging along any trajectory in the basin of attrac- 
tion. As mentioned before, entropy optimization would not converge for such a 
common invariant measure. 

To overcome the obstacle of singular measures, we propose a three step 
process: (?) regularization, (ii) entropy optimization, and (iii) inversion. Reg- 
ularization conditions the moment sequence into a feasible input to the en- 
tropy optimization, converting the original moment sequence into moments of 
a bounded, integrable phase function. The entropy optimization step can then 
be used to recover a closed expression for the phase function approximation. 
In the inversion step, point-wise evaluations of the phase function are used to 
recover an approximant of the original measure. 

The proposed regularization of the moment sequence of a singular, positive 
measure derives from an original idea of A. A. Markov to study the moment 



sequence through its complex generating function. We start by a simple observa- 
tion that an analytic function mapping a domain into the open upper-half plane 
admits an analytic logarithm whose imaginary part (the phase) is bounded from 
below by and from above by tt. The passage from a positive measure to the 
phase function through a canonical integral transform, obeying the above prin- 
ciple, has circulated in the Russian literature in connection with the century 
old works devoted to the one dimensional L-problem of moments. The early 
articles by M. G. Krein, N. Akhiezer and A. Nudelman on the subject offer a 
comprehensive account of this method [14, 15]. 

In the present article we go beyond one dimension, considering Fantappie 
transforms of positive measures supported by a wedge in M'^ [16, 17]. The 
existing methods of harmonic analysis on tube domains enter naturally into 
the picture offering to the maximum entropy reconstruction method a solid 
background. The much nicer sequence of moments of the phase function are 
obtained from the moment sequence of the original measure via a non-linear 
recurrent operation. A thorough investigation of the multivariate moment via 
asymptotic expansions of the Fantappie transform of the underlying measure 
was undertaken by Henkin and Shananin [18, 19], whose work we take as a 
basis for ours. 

While the entropy optimization provides a standard reconstruction proce- 
dure for the phase function, the approaches to inversion for one- and multi- 
variate problems are different. In one-dimensional case, we can make use of the 
well known Plemelj-Sokhotski formulas [20, 21] to complete the inversion step. 
The formulas, however, are difficult to generalize to multivariate settings [22]; 
instead, we propose a ray beam disintegration, based on a refined and partially 
forgotten one-dimensional analysis of the phase regularization due to Aronszajn 
and Donoghue [23]. The ray beam approach reduces the problem to a setting 
similar to medical tomography, based on inverse Radon or Laplace transform 
methods [24. 25]. We believe this will be a fruitful approach that we plan on 
exploring in follow-up papers, so we only draft it in this paper. 

The paper is organized using the following outline. Section 2 briefly intro- 
duces the multivariate moment problem and the entropy optimization, including 
an example illustrating lack of convergence for a Dirac-(5 measure. In Section 
3 we expose the elementary aspects of the entropy optimization method, in the 
case of one real variable for unbounded and bounded supports, using, respec- 
tively, power moments and trigonometric moments, i.e., Fourier coefficients. 
Section 5 is devoted to generalization to multivariate problems, through the 
phase regularization of the Fantappie transform of a measure supported by a 
wedge in Euclidean space (Section 5.1) and by special compact domains in Eu- 
clidean space (Section 5.2). A Riesz-Hcrglotz formula is derived, in the spirit of 
[26, 27], with a couple of examples on product domains. 

The present article remains at a theoretical level, leaving for a continuation 
of it to deal with further practical aspects: numerical experiments, the error 
analysis and examples from dynamical systems. We do, however, present prac- 
tical algorithms that are essential for moment conditioning, the Miller-Nakos 
algorithm in Appendix A, described in [28], and a recent algorithm for entropy 



optimization, described in [4], in Appendix B. 

We dedicate this work to the late Israel I. Gohberg, legendary figure of mod- 
ern operator theory and function theory. His original and highly influential ideas 
have permanently shaped moment problems and the entropy method referred 
to in the following pages. 

2. Preliminaries 

Let d > 1 be a fixed dimension and let K be a closed subset of the Euclidean 
space M''. Fix a finite set ^4 C N'' of multi-indices. The truncated moment 
problem with supports on K and monomials labeled by A consists in finding (as 
effectively as possible) a positive measure fi supported by K, with prescribed 
moments 

7a = / x^dnix), a e A. (1) 

Jk 

In case the set K is unbounded, it is implicit that the above integrals converge 
in Lebesgue sense. Throughout this article we adopt the multi-index notation 

x" =x'^'x'^\..x'^-', xeM.'^. 

A few basic questions are in order: 

1. Characterize all sequences of moments {act)aeA associated to positive mea- 
sures carried by the set K . 

This question can be rephrased in terms of the formal integration functional 

aGA a€A 

Let us denote by M[.T]yi the linear span, in the ring of polynomials R[x], of all 
monomials x" , a € A. 

A necessary and sufficient condition that a linear functional L : R[x] — > M. 
is representable by a positive measure supported by the set K is that L in non- 
negative on all elements / G M[x] which arc no n- negative on K . Then L can be 
extended via a Hahn-Banach construction to a positive linear functional on the 
space of continuous functions on K, with polynomial growth at infinity. This 
observation remains however of a limited theoretical importance, and it becomes 
effective only when simple characterizations of non-negative polynomials on K 
is available. Fortunately, in the case when K is a basic semi- algebraic set, such 
" Positivstellensatze" were recently resurrected and a good collection of examples 
is available, see [29]. 

The single variable case is the simplest and best understood. The following 
result goes back to Marcel Ricsz [30]. 



Theorem 1. Let n be a fixed degree and {a,h) an interval on the real line, 
hounded or not. A positive measure fj. carried by the closure of (a, b) exists, with 
moments 

Ik = / x''dfi{x), < k <n, 



> / x"dfi 



if and only if the associated functional L satisfies L{f) > for all polynomials 
f{x) ~ Cq + cix + ... + CnX^^ which are non-negative on (a, b). 

Three cases are distinguished, and they correspond to classical moment 
problem studies: (a, 6) = (0,1), known as the Hausdorff moment problem, 
(a, b) = (0, oo) known as Stieltjes moment problem, and (a, b) — (— oo, oo) 
known as the Hamburger moment problem. In each separate situation a full 
characterization of all non-negative polynomials on (a, 6) is available, with the 
result of making the above M. Riesz result effective. We refer the reader to 
[31, 32] for full details. 

2. Knowing that problem (1) is solvable, find constructively one particular 
solution. 

As a general rule, any attempt to solve the truncated problem (1) starts 
with the observation that the set of all solutions 



{^ > 0; / x"d^i = 7a, a e A] 

JK 



IK 

is convex and closed in the weak-* topology. If we include a = among the 
elements of the index set A, then all elements of E have fixed total variation. 
Thus, in this case, on a compact support K , the set of solutions E is compact 
in the weak-* topology of all measures. 

Among all elements of the solution set E the extremal ones are the first 
to be detected by linear optimization methods. For example, in the case of 
the three classical truncated moment problems on the line, they correspond to 
convex combinations of point masses. Their support is identified with the zero 
set of orthogonal polynomials, and the multipliers of the Dirac measures are 
also computable in terms of the diagonal Fade approximation of the series: 

_70 71 In 

z 



y2 '" -,n+l' 



Stieltjes original memoir remains unsurpassed for a careful analysis of this ap- 
proximation scheme, see for instance [31]. A basic observation in this direc- 
tion, providing an extremal solution to Stieltjes moment problem with the data 



(707 ■■■il2n-i) is the following: assuming that the Hankcl matrices 
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71 72 • ■ • 7n 



\7n-l 7n 



/71 72 

72 73 



72n-l/ 



\7« 7 



n+1 



7n \ 

7n.+l 

72n-lJ 



(2) 



are positive definite, the one step completion (70, •■•,72n-i,72Ti) so that the 
determinant 

70 71 •■• In 

71 72 • ■ • 7«+i 
. =0 

7n 7n+l • ■ • 72n 

vanishes, has a unique, necessarily finite, atomic solution. 

Since the computation of the roots of an orthogonal polynomial is not 
friendly from the numerical point of view, the search for other special solu- 
tions of the truncated moment problem led to adopt a statistical point of view, 
and consider "the most probable" solutions, with respect to a non-linear, con- 
cave functional. Recent applications (in particular to continuum mechanics) use 
to this aim the Boltzmann-Shannon entropy, see [1, 3, 7, 33-35]. The entropy 
maximization method for the trigonometric moment problem stands aside for 
clarity and depth in this framework, see [36]. 

A great deal of recent work, cf. [3, 9], has clarified the existence of maximum- 
entropy solutions, especially in some degenerate cases. We start from there, and 
add a computational/numerical analysis component to the study. 



3. Maximal entropy solutions in ID 

For the sake of clarity we digress and specialize the above discussion to the 
simplest and best-understood framework. Namely, we discuss below the exis- 
tence and uniqueness of maximum entropy solutions to the truncated moment 
problem in the case of a single variable. 

3.1. Basic properties 

Although an abstract, fairly general treatment of the maximum entropy 
method is nowadays available, see or instance [3, 34], we specialize below on 
an interval of the real line. To this aim, we go back to M. Riesz' existence 
theorem stated in the previous section. Namely, n is a fixed degree and (a, b) 
is an interval on the real line, bounded or not. We start with the moment data 
7o, ..., 7„, and seek a positive measure /x carried by the closure of (a, b) satisfying 



7fc 



x''dfi{x), < k < n, 



and 



> 



x"dfj.. 



We search fj, of the from d^(x) = exp(Ao + \ix + ... + \nX^)dx, assuming 
that the integrabiUty condition 

b 

exp(Ao + Aia; + ... + Xnx")dx < oo 

is assured by the choice of the parity and sign of the leading term. For instance, 
in case a = 0, 6 = oo we must have Ap < and Ap+i = Ap+2 = Ki = 0; or in the 
case a ~ — oo, 6 = oo we must have A2p < and A2p+i = A2p+2 = A„ = 0. We 
denote by A (by omitting the subscript n) the set of all such multipliers which 
produce integrable exponentials. 

The proper choice of the parameters Afc is made by imposing the optimality 
(maximum entropy) condition: 

sup<^ Ao7o + ... + A„7„ - / cxp(Ao + Ai.T + ... + A„a;")dx ^ (3) 

where the supremum is taken over all admissible (i.e. integrable exponen- 
tial) tuples A = (Ao,...,A„). Let us similarly denote 7 = (70,..., 7„) and x = 
(l,x,x^, ...,a;"), where the latter is considered as a variable point on the Vero- 
nese curve described by the list of the first monomials. 

The starting point of our discussion is the observation that the functional 

L : A — > R, i(A) = A • 7 - / exp[A • x]da;, 

J a 

is concave. Indeed, whenever the partial derivatives are defined (for instance in 
the Euclidean interior of A), we have 



dXidXj 



f'' ■ 

/ a;*"'""' cxp[A • x]dx. 

J a 



In the above Hessian, we recognize the negative of the Hankel matrix of a non- 
atomic positive measure, whence the strict negative definiteness. Moreover, 
the inner critical points of the functional are given by the vanishing gradient 
conditions: 

dL 



^r^' 



/ x'-' exp[A • xjda; ~ 0. 

J a 



The difficulty related to the described method lies in the complicated struc- 
ture of the set A of admissible multipliers. While for a bounded interval (a, h) 
this set is the whole Euclidean space A = ]R"+^, the case (a, h) = (0, 00) requires: 

A = [M" X (-00, 0)] U [M"~^ X (-00, 0) X {0}] U ... U [M x {0} x ... x {0}]. 

And similarly when (a, b) = (—00, 00). On the positive side, we remark following 
Junk [3] that in all cases the assumption that 7 is a moment sequence implies 

lim L{\) = —00. 



Thus, in the bounded interval case, the optimization problem (3) always 
has a solution, and by strict convexity, this is unique. Note that in this situa- 
tion, the positivity conditions in M. Riesz Theorem (or equivalcntly HausdorfF 
finite difference conditions) are necessary and sufficient for the existence of an 
exponential type solution to the truncated moment problem, see also [1] for a 
detailed discussion. 

A much more delicate analysis is required in the case of Sticltjes moment 
problem (a, b) = (0, oo). For this case it is very possible that the extremal value 
in problem (3) is attained on the boundary of the set A. Assume for instance 
that 

sup < Ao7o + ... + A„7,i - / cxp(Ao + Xix + ... + Xnx")dx > = 



/ cxp((To + Aix 

J a 



a-o7o + •■• + o-„7„ - / cxp(cro + Aix + ... + cjnx"-)dx 

J a 

where a = (ctq, ..., ct„) G A \ int A. That is, there exists an index Q < p < n with 
the property 

CTp-i < = CTp = ... = cr„ 

if p > 1, or simply 

= 0-1 = ...= (Jn 

in case p = 1. Anyway, then only lateral partial derivatives ^^(c) exist for all 
P ^ j ^ n. Since cr is a global maximum, we infer 

J 7j - /o°° x^ exp[CT • x]dx = ■|^(cr) > 0, p<j<n, 
\ l3 - /o°° ^^ exp[cr • x]da; = |^ (cr) = 0, j < p. 

Note that above, the exponential density depends only on p parameters 
(tToj ■•■:0'p-i)j whence it is normal to expect that only the first p moments are 
matched. 

A detailed analysis of the decision tree resulting from the above observations 
goes as back as 1977 to Einbu [37] and it was much clarified in the recent works 
by Junk [3] and Hauck, Levermore and Tits [9]. We reproduce below, following 
Einbu and Junk, the main phenomenon, in the form of an analysis of a one step 
extension. 

Suppose that, for the truncated version of Stieltjes moment problem, the 
initial segment of moments 

(7o,7i,---,7"-i) 

is realized by the maximum entropy method, that is there is an admissible tuple 
a = ((To, ••■, cTn); such that 



7j 



/>oo 

/ X-' exp[cr • xjcfx, < J < n — 1. 

JO 



This implies that Hankel's positivity conditions (2) hold true, and that the 
(lateral) partial derivatives of the function L(A) vanish at A = cr. 

We assume next that the extended moment sequence (70,71, ...,7„_i,(5) is 
also realizable by the maximal entropy method. Hankel's positivity conditions 
(2) imply 

<^ > 7n(niin), 

where the bound 7„(min) is a rational function of the data (70,71, ...,7„_i), 
expressed as a quotient of Hankel type determinants. Define 



7„(max) 



/>oo 

/ x" exp[(7 • x\dx. 
Jo 



This corresponds to the boundary point (cto, ..., ct„, 0) S A„, and in addition we 
know that the function L : A„ — !■ R, when restricted to A„_i x {0}, has null 
partial (lateral) derivatives at (ctq, ...,cr„,0). Assume that 



7j 



/"OO 

/ X-' cxp[r ■ x]dx, < j < n, 
Jo 



where r G A„. In particular r„ < 0, or t„ = 0, in which case, by the uniqueness 
of the maximum entropy solution t = (ctq, ..., cr,i, 0) and S = 7„(max). 
Assume that t„ < 0, so that 



7j - / x' exp[T • x]dx = ;^t-(t) = 0, 
Jo c'-^j 



where jn = S- Thus r is a global maximum for the function L defined on A„, 
and in particular L{t) > L{(7o, ...,an,0). By analyzing the restriction of the 
concave function L to the linear segment joining inside the set A„ the points r 



and ((Toj •■•>o'n, 0) wc infer g^ 



3L(tr+(l-t)(go,...,g-„,0)) 



< 0, or in other terms 



6 < A„(max). 



In conclusion, assuming that the finite moment sequence (70,71, ■•■,7n-i) is 
representable by a maximum entropy solution of the same degree, the extension 
(70,71, ...,7„-i,7„) has the same property only if 

7„(min) < 7„ < 7„(max). 

One step further, when investigating only the solvability of Stieltjes problem 
with data (7o,...,7n) by the maximum entropy solution without assumptions 
on the projected string (70, ...,7„_i), the upper bound 7ji(max) may become 
infinite, see for details [3]. 

3.2. Recurrence relation for the moments of an exponential weight 

The maximum entropy method for solving the truncated moment problem 
invites us to have a closer look at the full string of moments of an exponential 



of a polynomial weight. We enter below into the details of these computations, 
in the case of Stieltjes moment problem. 

Fix an integer n > and consider the polynomial 

P{x) — (To + aix + ... + (T„x", 

with real coefficients and cr„ < 0. Denote by 

/•CJO 

7fc = / x'' exp[P{x)]dx, fc > 0, 
Jo 



the moments of the density e^^'dx. An integration by parts yields, for all 
A;>0: 

/■oo ^fc+1 /-oo fe+1 /-oo k+l 

Hence, a finite difference equation relates every string of n + 1 consecutive 
moments: 

(fc + l)7fc + cri7fc+i + 2cr27fc+2 + ■ • ■ + nCT„7fe+„ =0, fc > 0. (4) 

Since ct„ 7^ 0, we obtain the following simple observation. 

Lemma 2. Let P{x) he a polynomial of degree n, with negative leading term. 
The moments of the density e^^^'dx are recurrently determined by (4) from the 
first n moments. 

Specifically, the linear dependence 

fc + 1 CTi {n ~ l)o-„-i 

7fe+n = 7fc 7fe+l - ■ ■ ■ 7fe+n-li 

nan na„ nan 

holds. By changing the running index, we find for all m > n : 

m — n + I a\ C*^ ^ l)c,i-i 

7"i 7Tn — n 7^?i— n + l ■ ■ • 7Tn— 1- 

nan nan nan 

Let M' = max"-^ ^^^ and M = max(Af', ^i^ ), so that 

max|7j| < (- + nAf ) max |7j| . 

j<ra \nan\ j<rn—l 

Therefore there is a positive constant C and a positive integer A^, such that 
max |7j| < C""(to + A^)!, 777 > 0. 

j<m 

Consequently, Stirling's formula implies 

lnmaxj<,„|7j| (777 + A^)(ln(m + A^) - 1) ln(27r(777 + TV)) 
< L' H \ ^ , 

777 777 2m 

10 



and in particular 

In Tna"X"„-^™ ln/„ I 



ln.maxj<m\jj\ ^ ^, , ln(m + A^) 



m m 

where C" is a positive constant. 

In conclusion, there is a positive constant 7, such that 

00 _. 00 _. 00 -. 

Z^ I |i/m - Z^ [niax,<m hiW^'"- -^ ^ m + N ~°°' 

m=0 17™ I m=0 ^ •'-™ ' '-"J m=0 

According to Carleman's uniqueness criterion (sec for instance [31]) we obtain 
the following result. 

Theorem 3. Let P{x) be a non- constant polynomial with negative leading term. 
Then the moment problem with density e^dx is determined. 

We translate this statement for the reader who is not familiar with the 
terminology: if a positive measure fj, on [0, 00) has the same moments as e^dx, 
then /i = e^dx. 

3. 3. Existence 

We have seen in the previous sections that not every truncated sequence 
of moments (70, 71, ..., 7„) on the semi- axis can be achieved by the maximum 
entropy method, within the same degree. That is, it is not true that there 
always exists an admissible polynomial P{x) of degree n or less, such that 

7fe = / x'^e^^'^Ux, k<n. (5) 

/o 

To give the simplest example, consider the sequence 

70 = 1, 71 = 72 = • ■ ■ = 7n == 0- 

Obviously, the Dirac mass 5q has these very moments. However, there is no 
polynomial P, of any degree, such that 

/•OO 

= 71 = / xe^^'^Ux. 
Jo 

Simply because the integrand is non-negative and non-null on the interval of 
integration. 

Our study is motivated by the need to solve this pathology. In the follow- 
ing sections we indicate a method to overcame the limitation of the maximum 
entropy method to absolutely continuous measures. Along the same lines, some 
recent works proposed different regularizations, see for instance [38] 
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4. Single variable: Conditioning using the Cauchy transform 

The recent works of Junk [3] and Hauk, Levermore and Tits [9] clarified 
which positive densities p are appropriate for the maximum entropy reconstruc- 
tion method. A thorough analysis of the convex structure of the truncated 
moment set of these distributions, e.g., extreme points, facets, was carried out 
in the cited works, with significant applications for the kinetic theory of gases. 
In particular, singular measures are especially poor candidates for maximum 
entropy reconstruction, as seen from the example in Section 3.3. In an attempt 
to enlarge the class of measures for which such well established reconstruction 
methods work, we propose a regularization procedure which will produce an 
admissible input for the entropy optimization procedure. 

The goal of our procedure is to reconstruct a possibly singular measure /u 
by transforming it to a continuous measure (j){t)dt, whose density (f> we term the 
phase function. The entire measure reconstruction procedure can broken down 
into three steps: 

1. regularization based on moment data of ^, 

2. density reconstruction (using entropy optimization) of cf), 

3. inversion, i.e., recovering a measure /x* « /i, from point-wise knowledge of 



We stress here that it is not our aim to improve on the density reconstruction 
procedure, i.e., the entropy optimization, itself. Rather we focus on moving the 
density reconstruction where it can be performed with assured convergence, by 
inserting the regularization and inversion steps. It could be very well possible 
that other density reconstruction methods, e.g., basis pursuit, wavelet-based 
reconstruction, could be used instead of the maximum entropy for the general 
reconstruction problem, however, we do not explore these options here. 

In this section, we first focus on measures whose support lies in a one- 
dimensional space. In this case, the entire procedure is based on a simple idea 
of A. A. Markov [39], widely used in function theory, employing Cauchy trans- 
forms. Cauchy transforms serve as an analytic tool to study complex generating 
functions of the moment sequences. The regularization step is based on repre- 
sentation theorems for the generating function of the moment sequence, while 
the inversion step is grounded in Plemelj-Sokhotski formulas, which can be used 
to reconstruct the densities on the original domain. When the domain is one- 
dimensional, Plemelj-Sokhotski formulas can be formulated through a Hilbert 
transform, which is easily evaluated numerically. Therefore, such a reconstruc- 
tion results in an algorithm that can easily be implemented in a computer code. 

In Section 4.1, we first give the procedure for measures with arbitrary sup- 
ports in K, based on power, i.e., monomial, moments of the measure ^ as input 
data. If the support of measure is contained in a compact interval, we can em- 
ploy trigonometric moments instead, which are preferred numerically to power 
moments. The regularization procedure for compact supports is developed in 
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Section 4.2, and is somewhat more technical than for the unbounded case, yet 
the spirit is the same. Based on insights for one-dimensional domains, in Section 
5 we discuss how the procedure might be extended to measures supported in 



4.I. Unbounded support 

Define the Cauchy transform of a measure fi, with support in M, as 

Ci^) = / ^. (6) 



Markov's observation is the following: assuming all integrals exist, the Cauchy 
transform of a positive measure on the line is of Nevanlinna class, i.e., it has a 
positive imaginary part in the upper-half plane: 

CKz) ~ Cfi{z) _ f Qzdfi{x) ^ p c^2 > 



2* Jb |a; — ^1 

Hence the phase Qhi[C^{z)] is a harmonic function in the upper half-plane, 
uniformly bounded from below by zero and from above by tt. The boundary 
values along the real line of 51n[C/i(z)] produce an integrable, positive and 
bounded density 0, satisfying: 

1 + C^i{z) = cxp / ^^^, 32 > (7) 

Jr X- z 

i.e. 

1 + Ci^l{z) = cxpC(j){z), 3z > 0, (8) 

where we slightly abuse the notation when we use C4>. The dictionary between 
properties of fi and density cj) was established by Aronszajn and Donoghue [23]. 
The most important, of course, is the existence and boundcdncss of (j). As is 
bounded even if ^ is singular, we consider 4>{t)dt to be a regularization of dfi. 

Practical benefit of this expression comes from the ability to use it without 
knowing the closed-form expressions for measures involved. The Cauchy trans- 
form is the (complex) generating function for moments of /i, i.e., its expansion 
at z = oo is given by 



(C/i)(^) 



n=0 ^ 



where a^(n) — j^f''diJ,{t), and a^{n) defined analogously. ^ Solving for C4> and 
using the series expansion ln(l + z) = — X]i^i(~l)"^"/'^ yields the following 



■^Thc non-linear transform of the moment sequence was exploited in the theory of the phase 
shift of perturbed spectra in quantum mechanics, see [40, 41]. 
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equality between power series: 

oo 

C<t>{z) = ln[l + C^l{z)] = - ^ T[-C^^{z)]' 



fe=i 



oo / N oo _, 

E acf,(n) _ v^ 1_ 

n=0 A:=l 



.n=0 



The Miller-Nakos Theorem [28], whose complete proof we bring in the Appendix A. 
gives a recursion for evaluation of moments a^{n) from moments a^{k) for 
fc = 0, . . . ,n, 

where [5jv(2)]|r indicates the coefficient next to z~^^^^\ in the fc-th power 
of the truncation Sn{z) = X]n=o Q/j('t-)-z~^"^^^ of the generating power series. 
Such a triangular property is essential for practical problems: we will typically 
have access only to truncated moment data and we do not wish to establish any 
a priori ansatz, especially not a^{n) — for n > N. 

At this point, we have set up moment data such that most density recon- 
struction procedures apply: density (j> is bounded and compactly supported, 
making it possible to reconstruct it using entropy optimization described in 
Section 3. Such a procedure produces an approximant 

N 

(j)* {x) = exp ^ akx'' (10) 

k=0 

that converges to density cf) as the number of available moments A^ increases. 

The inversion step describes how the point- wise knowledge of approximant 
(^* « (/) is used to compute an absolutely continuous measure /i* that approxi- 
mates the original measure fi. In this paper we do not claim to obtain quanti- 
tative convergence results on ^* — > /i, especially when /i is a singular measure, 
however, we stress that, for singular measures, a classical reconstruction pro- 
cedure like entropy optimization might not produce any results. Therefore, we 
view our results in this paper as a starting point for further investigations of 
approximation of singular measures. 

To a smooth entropy optimizer (j)* corresponds a measure n* with a density 
p — dfi* /dx, which we use to approximate the original measure fi. The lynchpin 
of the inversion procedure, i.e., evaluation of p from knowledge of (p* , is the 
existence of boundary limits linie_>o C/u(x ± ie), for e > 0. The limits exist 
independently pointwise, and, assuming that p G iy^(R) is of Holder-class, the 
Plemelj-Sokhotski formulas, e.g., [20, §14.11] or [21, §3.7], establish that it is 
possible to evaluate p pointwise from Cauchy transforms of p as 

p(x) = - — lim[C/i*(.T + ie) — Cp*(x — ie)]. 
2m eio 
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As limits exist independently, and exp is analytic, we can formulate them in 
terms of analogous limits for Cauchy transforms of the phase function Ccfi* , i.e. 

by (8), 

p(x) ~ - — limfcxp C(/)* (a; + le) — exp C(j}*{x — ie)] 



1 

27ri 



exp lim C0* (x + ie) — exp lira C(j)*{x — ie) 



Moreover, the Plemelj-Sokhotski formulas provide explicit expressions for each 
limit: 

lim ^Ccj,*{x ± le) = ±]-cj,*{:x) + ^H0*(x), (11) 



where the Hilbcrt transform is 



-KJ t ~ X 

It follows then that the expression for p is given by: 

p{x) = — exp [— 7r?^0*(x)] sin7r(^*(a;). (13) 

TT 

This formula connects density p = dp/dx with the phase density function 0, or, 
in the case of moment closure by entropy optimization, a smooth approximant 
(f)* to the phase density function 0. 

The formula (13) is numerically practical: the entropy optimization provides 
us with a closed formula for </>* , while its Hilbert transform is easily numerically 
evaluated via the Fast Fourier Transform algorithm. Therefore, here we have 
obtained a practical inversion formula for approximating a singular measure // 
via an absolutely continuous measure p* with density p = dp* / dx. 

4-2. Compact support 

When the measure p is supported on a known compact interval, we can use 
trigonometric moments, instead of power moments, in the process given above. 
The resulting process is more numerically robust, as trigonometric functions 
are orthonormal and bounded as a family, unlike the family of monomials on an 
arbitrary interval. 

Let /x be a measure on the interval A = [— tt, tt) that induces the measure 
p, on the boundary 9D of the unit disk B) C C A known relation is then 
dp{9) = —i(dp{(), for ^ = e*^ S dD. Define the circular Cauchy transformation 



ICpiz)^ ' f '^^(^^ 



27r Job Q- z 

for z ^ (9D. Using the equivalent arc-length formulation clarifies the difference 
between 

Kpiz) - ^ r MO) 



2tt J_^ 1 - e-'" 
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and the Cauchy transform on the Ime C/x(z) = J^^ diJ,{x)/{x — z), cf. (6). 

The function K,^ is holomorphic inside intD, /C/i S 0(D), where it has the 
Taylor expansion 

CJO 

JCfi{z)=^Y,r^,{k)z'', 

k=0 

with complex trigonometric moments 

r,(fc)--^/c-^ (14) 

^TT Jan < 

serving as coefficients. 

The imaginary part 5/C/i(z), is positive for positive measures, as the imagi- 
nary part of the kernel is 

1 / z -i \ 1- ^(e-'^z) 



Ani \1 - e-'" z l-e^^zj 27r |1 - e-*^z|^ ' 

and ze"*^ < 1 when z G intD. Consequently, the argument oi ICfi{z), with the 
appropriately chosen branch of the logarithm, 

F{z) ^ -i InlCfiiz) e 0(B), (15) 

is of Caratheodory class: it is a positive function, with a bounded real part 
5RF(z) £ [0,7r], which corresponds to the bounded angle of /C/i(z). 

The following classical theorem allows us to obtain a representation of Caratheodory 
class functions in terms of bounded densities on a circle [e.g. 42, §12.10]: 

Theorem 4 (Riesz-Herglotz). Let F e 0(ID)) be such that ^F{z) G [0,c] for 
some fixed c > 0. Then there exists a function (p G L^{dl])) for which 

F{z)=i^F{0)+P(j){z), 
where 0(C) £ [0, c] pointwise and 

is the Poisson integral of (f){9) = (f){e ). 

The Poisson integral Vip can be rewritten in terms of the circular Cauchy 
transform K,4>: 

271" JdnC- z < 

where, again, T<^(fc) are trigonometric moments of 0, defined analogously to (14). 
The Riesz-Herglotz formula then reads 

F{z) = -i2/C(/)(z) - T0(O) -I- iQF{0). 
16 



We can compute the constants in the formula by evaluating it at z 
comparing it to evaluation of the definition (15) at the same point: 

F(0) - z3F(0) - T0(O) - %2{%T^(iS)\ - r^(0) + i3i^(0) 



and 



concluding that 



F(0) = -2lnzT^(0) = - -zlnr^(O), 



r40) = -, 3F(0) = -lnr^(0). 



Substituting these constants into the Riesz-Herglotz formula, and using the 
definition of F{z), we obtain the exponential representation of JC^{z): 



K:^l{z) = -ir^(O)exp[2/C0(z)]. 



(16) 



To compute the moments of (p, we relate the Taylor expansions of the func- 
tions above, and use T0(O) = 7r/2 to obtain 



^T^,{n)2 



exp 



oo 

2i^r^(n)z'^ 



where ff^(n) = r^(ri,)/Tp(0). As before, we use the expansion ln(l + z) 
— X]^i(~l)"'2;"/n to relate the series through expression 



k=l 



T^{k)2 



2^ k 



1 k 



.n=l 



T^{n)2 



A finite number M of trigonometric moments r^(fc) can then be computed 
using the Miller-Nakos algorithm (see Appendix A) when M moments T^(fc) 
are known. 

To invert the procedure, we assume that to approximate (f>, we are given 
a smooth density </>* : [— tt, tt] — >■ M, which corresponds to a continuous /x* 
with density p : [—tt, tt] — > R, i.e., diJ*{9) = p{9)d9. Density p can be evaluated 
point-wise using Plemclj-Sokhotski formulas (see Dynkin's chapter, section §6 in 
[43]), which evaluate non-tangential limits i-limj_i.2 /C/i*(^) and e-limj_^z /C/i*(^) 
at z e 9D. with the argument in domains ^ e int D and ^ G C/D, respectively. 
The A. Calderon's theorem asserts existence of such limits for cj)* S L^((9D). 

Due to analyticity of exp in (16), we can evaluate the non-tangential limits 
of ICp* , in terms of non-tangential limits ICcj)* of <f>*{C)d-C^ 



i-lim/C/i*(^) 
e-lim/C^*(C) 



-ir^(O)exp[2i-lim/C0*(O] 
-ir^ (0) exp[2 e-lim /C0* (^)] 



•^Wc slightly abuse the notation when we use K<f>* 
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Privalov's Lemma establishes that the non-tangential limits satisfy Plemelj- 
Sokhotski formulas for z G c'D: 

i-lim/C0*(O-Q0*(^) + ^0*(^), 

e-lim/Cr(C) = Q0*(z)-^0*(z), 
where Q indicates the singular integral 

0*(C)dC 



Q<\> 






C 



with analogous expressions holding for ICfi*{z), with density p{9) ~ dfi*/d9 
instead of cj)* . Therefore, to evaluate p{z) on 9B we seek the difference between 
the non-tangential limits: 



p{z) = -T^(0) <^ cxp[2i-lim/C0*(O] - exp[2e-lim/C0*(e)] 
= -i2T^(0)exp[2Q(/>*(z)]sin(/)*(z) 

The singular integral Q<f>*(z) can be evaluated on z = e*^ using the circular 
Hilbert transform of 0* : 

nA*f ^ 1 i ^ICK » f C + z riQdC , » ,^, 



/71 



IIT 



*^a)da + — 



= -£^cot^0(a)da + - 
where z = e'^. The circular Hilbert transform is, by one convention, 

Finally, substituting this expression into p{C,) = p{9), we get the evaluation of 
the density p{9) as 

pi9) = 2T^iO)cxp[n(^*i9)]sm^*{9) 

Practically, Hilbert transform is easily evaluated on a fixed grid using numerical 
Fourier Transform, e.g., FFT, so this formula can be employed when we have 
access to the evaluation of </>* on a fixed grid in [—it, tt]. 
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5. Several variables: Conditioning using complex Fantappie trans- 
forms 

The reminder of this paper deals with generalization of the regularization 
procedure to measures of several variables. To do so, we will replace Cauchy 
transform with a Fantappie transform, which is usually defined as a real inte- 
gral transform, relying on two sources: the harmonic analysis on tube domains 
over convex cones [44] and the complete monotonicity results, a la Bernstein, 
characterizing the Laplace and Fantappie transforms of positive measures on 
convex cones [18]. For expository material on Fantappie transform, see [16, §3]. 

We propose three different ways to complexify it, in order to use the general 
Riesz-Herglotz representation theory and obtain analogs to the phase function </>. 
The choice of the complexification procedure is based on a trade-off: presently 
we are able to obtain either theoretically general results with little practical 
value, or practically useful results which do not allow for as much theoretical 
breadth. We expect that the future research will bridge this gap between theory 
and computation. 

Take a solid, acute, closed convex cone F C M'', and its associated polar cone 

r* = {x(ER'^; uj-x>0, uj gT}. 

The cone F* will carry the support of the measure /x, while F will play the 
role of the "frequency parameter" , to use the language of signal processing and 
applied Fourier/Laplace analysis. The following characterization of the real- 
valued Fantappie transform is due to [18]: 

Theorem 5 (Henkin-Shananin). A function $ : (0, oo) xF -^ M. is the Fantappie 
transform 

$(-0,^)=/^^, (17) 

of a positive measure ^ supported by T* if and only if $ is (i) continuous, 
(ii) completely monotonic^ , and (Hi) homogeneous of degree —1, i.e., 

^{Xljq, \u}) = A~^<i>(a;o, w), wq > 0, a; G F, A > 0. 

To extend the Fantappie transform to complex domains, one has a choice of 
complexifying the offset parameter wq, the normal parameter uj, or both. We 
start with the full generality in Section 5.1, complexifying both loq and uj to tube 
domains, and develop the full regularization procedure, at a cost of providing 
no inversion formulas. Next, we constrain ujq = 1 in Section 5.2, obtaining a 
restricted tube domain, to provide some practical regularization formulas for 
measures on familiar compact domains, which stand in direct analogy to one- 
dimensional problem. A future research direction could explore Clifford algebras 



* A function <& is completely monotonic if it satisfies inequalities (— l)*^Djj ■■■^C.k ^(p) — ^' 
Vfc > 0,p G intr, where D^^, are partial derivatives along eoordinates ^i, ...,^fc £ F. 
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as a setting for generalization of the Plemelj-Sokhotski formulas, which were 
used to complete the inversion process in the one-dimensional case in Section 
4. A shorter, perhaps more immediately practical, procedure is given in Section 
5.3, where we treat w as a parameter, complexifying only cjq- As a consequence, 
we recover a single variable procedure at each value of w, which we term the 
partial Fantappie transform. The family of solutions, parametrized by w, could 
be used in a tomographic procedure to recover an approximant to the original 
measure ^. 

5.1. Unbounded supports and tube domains 

The Fantappie integral transform extends analytically to the complex do- 
main: 

$(uo, u)= f _M^^ sfty £ int r, Kuo > 0, 
Jy» uq-\-u- X 

retaining, by definition, homogeneity of degree —1 in the complex argument 
(uo, It) G C X C''. Let ri = (0, oo) x int F C W^^^ be the interior of the domain 
of continuity for the real Fantappie transform. 

The associated tube domain is the set Tji = E + ifi, where E = M'^'+^. Due 
to the different role played by the first axis, we denote elements by [zq^z) = 
((To -|- iojQ,o + iuj) S Tq. With this notation^ the domain of analyticity can be 
written as —iTn, i.e., $ G 0{—iTn). For clarity, we will use u to denote elements 
of —iTn, and z for elements of T^, with the obvious change of coordinates 
u = —iz = CO — ia, when z ^ a + iuo. 

Notice that when {uq,u) G —iT^, 

3fi(f>(wo,'«) ^ / — ^dfi{x) > 0, 

Jr* \uo +u- x\ 

following from the definition of the polar cone F*. Therefore, function i$(wo, u) 
is analytic, and has a positive imaginary part. It follows that its complex phase 
lni$(uo,'u) is well defined on —iT^, with the property 

31nz<i>(Mo,'") G (0,7r), {uq,u) e -iT^. 

Converting this expression to the tube domain, define F{zq,z) on the tube 
domain Tq by setting 

F{zt),z) = —ihii<^{—izQ,—iz), 
with a further simplification 

F(zo,^) = -iln[-$(zo,^)], (18) 



^Such a choice conforms with the existing conventions of harmonic analysis, and we apol- 
ogize in advance for all the resulting multiplicative imaginary unities. Fantappie transforms 
in the later sections simplify this convention. 
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due to homogeneity of $, or 

$(zo,z) = -expiF(zo,z). (19) 

Defined this way, the function F{zq, z) is analytic on Tft, and satisfies Ki^(zo, z) £ 

These properties make it possible to reveal the structure of functions F via a 
straightforward generalization of Riesz-Herglotz formula for analytic functions 
of positive real part (see Theorem 4), which is our next goal. 

Let i7 C K'' be an open, acute and solid convex cone, with associated tube 
domain Tq = R"^ + ifl. The Hardy space H^{Tfi) is defined as the space of 
analytic functions F : Tq — > C, such that 



\F\^ = sup / \F{a + iuj)f da < oo. 
wen jRd 



By a celebrated theorem of Palcy and Wiener, H^{Tq) is the space of Fourier- 
Laplace transforms of square integrable functions defined on the polar cone. 
The following result characterizes real Fourier-Laplace transforms [18]: 

Theorem 6 (Bernstein, Bochner, Gilbert). A function F : ft — > M. is the 
Laplace transform 



F{uj) = / e-'^-''d^i{x), 
Jn* 

of a positive measure fi supported by f2* if and only if F is continuous on fl and 
of class C°° and completely monotonic in the interior int J7. 

The extension from the cone il to tube domain Tq, for / G L^{il*,dx), is 
given by 

for z € Tq. Function F then belongs to H^{Tq) and the map / i— >• F is an 
isometric isomorphism between the two Hilbert spaces. For a proof and an 
overview of the theory of Hardy spaces on tube domains sec [44] . 

The reproducing kernel of the Hardy space, also known as Szego's kernel is 



S{z,w)^ — i-T / e''(^-'">^da;. 



for z^w € Tfj. Remark the homogeneity property: 

S{Xz,Xw) = X-'^S{z,w), A>0. 

The reproducing property has the following form: if F G H^{Tn), then the 
boundary limits, still denoted by F, satisfy 

F{a) = liiR F {a + iuj), 
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where the limit exists in L^fR'') and 



F{z)= I S{z,a)F{a)da. 



We focus next on functions F G A{Tq) which arc analytic in Tq and uni- 
formly bounded and continuous on its closure. Then the function z t-^ S{z, w)F(z) 
belongs to H'^{Tq), for every w GTq, and 



S{z,w)F{z) 
and by complex conjugation 



S{z, a)S{a, w)F{a)da, 



Siz,w)F{w)= Siz,a)Sia,w)F{a)da. 
By adding the two identities we obtain 

^^^^^^F(z)+£H ^ f s{z,a)S{a,w)^F{a)d<j. 
The restriction to the diagonal of the above formula yields 
3fiF(z)= f P{z,a)^F{a)da, 
where 



P{z,a) 



\S{z,a)\' 



S{z, z) 

is Poisson's kernel. Again, see [44] for full details. 
Fix a point a G zil, and subtract the identities 



z 6 Ta, ae M^ 



S{z,a)[F{z) + F{a)] 



S{z,a) 



F{a) + F{a) 



2S{z,a)S{a,a)^F{a)da, 

i 

^(cr, a)| S{z,a) 



S{a, a) 



-^F{a)da. 



One finds 

S{z,a)[F{z)-i'iF{a)] = 



2S{z,a)S{a,a) — 



\S{a,a)\ S{z,a) 
S{a,a) 



^F{a)da, 



a tube domain analogue of the Schwarz formula, which relates the values of an 
analytic function (in the disk) to the boundary values of its real part. We call, 
by way of natural analogy to the similar integral kernel for the disk, 



H{z,w; a) — 2 



S{z,w)S{w,a) \S{w,a)\'' 



S{z,a) 



S{a, a) 
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the Herglotz kernel associated to the tube domain Ta, with z,w G Tq, a <E ifl. 

In particular cases one can obtain from here an integral representation of 
all analytic functions in Tq possessing positive real part, what it is customarily 
called the Riesz-Herglotz formula, see for details [27]. Fortunately, our aim is 
more modest, having to deal only with the Fantappie transforms of measures 
appearing in the previous section. 

From now on, we return to the convex cone $7 = (0,cx)) x intF C ]R''+^ 
appearing in definition of the real Fantappie transform (17). Specifically, the 
analytic function 

F{zo,z) = -iln[-$(zo,z)], (zo,z) e Tn = M^'+i + ^17 

satisfies 

0<3?i^(zo,z)<7r, {zo,z)eTn 

and, due to the homogeneity of $ it growths logarithmically along rays contained 
in Tn: 

|F(Azo,Az)| <|F(zo,z)| + |lnA|, (zq, z) e Tn, A > 0. 

Fix a point a G iQ and consider the translated functions 

F^{zo,z) = F[{zo,z) + ea], 

so that they are analytic on the closure of Tji and of logarithmic growth along 
rays. Due to the homogeneity of degree — (d + 1) of the reproducing kernel S of 
To, we deduce that for every e > 0, the function F^S & H'^{Tn), and in view of 
the computations above: 

F,{C)^i^F,{a)+ [ HiC,a;a)^Fei<j)da, C = i^o, z) e Tn, SiC,a) ^ 0. 

Note that ^F^ £ (0, n) on Tq, as restriction of the original function to a subset of 
the tube domain. By passing with e to zero and to a weak-* limit in T°°(]R''+^) 
we obtain a function e L°°(R'^+^), < 4> < tt, a.e., representing F as follows: 



T(C) = i5T(a) + / i/(C,a;«)0(<T)da, C-(^o,z) 



GTo, SiC,a)^0. 



In all instances of interest Szego's kernel S'(C, a) does not vanish at all points 
C, a G Tn, producing a genuine integral representation of F. 
We collect the above remarks into a formal statement. 

Proposition 7. Let T C W^ be a closed, solid and acute convex cone and let n 
be a finite mass positive measure supported by the polar cone T* . The Fantappie 
transform of the measure jj, admits the exponential representation: 

f ^^^^ = -cxp[zT(zo,^)], iz,,z)&Tn, 
Jr* zo + z ■ X 
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where Q = (0, oo) x intF. In its turn, the analytic function F admits the integral 
representation 



FiO ^iC+ HiC, a; a)0(a)da, C = (^o, z) G To, 

JK'i + i 

where (f) : M'^'+^ — > [0, tt] is a measurable function and C G M is a real constant. 

Wc assume in the above statement that a G iil is fixed and S'(C, a) ^ 0. 
In this sense the, possibly singular, measure /x is regularized by the absolutely 
continuous measure (j){a)da. The integral kernel ^(C, c; a) is of little practical 
use in its full generality. However, in particular cases, to be discussed in the 
rest of the article, the preceding Markov exponential transform regularization 
becomes more accessible. 

One question which naturally arises in the above statement is: is it possible 
to characterize the bounded densities (f> appearing in the integral representation 
of an analytic function F G 0(To) that satisfies < 5RF < 1? The answer is 
yes, but the conditions imposed on cj) are not friendly. They were discovered 
a long time ago, in the case of the polydisk [45] and more general symmetric 
domains [27] . We merely indicate these conditions in the case of non- vanishing 
Szego kernel and sketch the proof. 

Proposition 8. Let fl (Z M."^ be an open, acute, solid cone, let a G ifi and 
assume S{C,,^) 7^ 0, Ci^ ^ ^fi- ^'^ element (j) G L°°{M.'^) is the phase of an 
analytic function F eOiTn), < KT < 1 : 



F{C)=iC+ f H{C,u-a)cj>{u)du, C G Tn, (20) 

7 G M, if and only if < 4> < 1, a.e., and the ''moment conditions" 

JR'' 



hold. 



Proof. In order to prove the non-trivial implication, let </> G L°°{M.'^),0 < (p < 
l,a.e., and define the function F{Q by formula (20). In view of the definition 
of Hergiotz' kernel, by taking ^ = a we find 

Fia)-^C=n^^ct.iu)du. 
J b{a,a) 

whence by addition 

FiO^FW) = 2J'-^^^p^<t>iu)du. 
J S{C,a) 



Write this formula for F(^) + F(a), too. According to the moment conditions 
in the statement, we find again by addition 

F(C)+F(O=2|^i^^|Mj^0(.)du, C^eTn. 
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In particular, ^F is the Poisson's transform of </>, and therefore < 5RF < 
1. D 

5.2. Supports in special sets and restricted tube domains 

The reader should be puzzled by now by the way too abstract and useless 
level of this article. It is time perhaps for some examples of phase regularity 
applied to measures supported by three basic convex shapes in euclidean space: 
the orthant, the euclidean /2-ball and the ^i-ball. The first one will especially 
be dear to control theorists, because it contains in the particular case of one di- 
mension familiar computations of Laplace transforms. We include the euclidean 
ball as the commonly occurring domain for measures, and the h ball as it results 
in trigonometric moment data for phase functions, which is an appealing set up 
for entropy optimization. 

To remove some normalizations that encumber the computations we break 
with the generic convention Tfi = W^ + ifl. Instead, at the beginning of each 
example, we specify the domain, redefine the Fantappic transforms and then 
summarize the derivation which was detailed in the Section 5 to obtain the final 
result. 

5.2.1. The orthant 

Let fi = (0,00)'' be the open positive orthant in M'', self-dual in the sense 
fi* = f2, the closure of itself. Szegb's kernel of the tube domain over Vt is 

(271-)" Jn- {2Tny- ^J-^ Wk - Zk 

With the selection of the reference point a = {i,i, ...,i) ^ iil, Herglotz kernel 
becomes 

S{z,u)S{u,a) \S{u,a)\'^ 
H [z, u;a) — 2 



S{z,a) S{a,a) 

d 

2TT ^~'^fe TT- 

J--^ (uk ~ Zk)(uk + i) .J--'- 1 



(27r)'i [ J-^A (life - Zk){uk + i) ^J-^ 1 -t- ul 

2n-Lf_^ L_-)A_Lf_t 

J- J- 9777 \ 111. — y.i^ II 1, -\- i I -^ J- O.iri \ in. - 
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2TTi\Uk~Zk Uk + iJ f-\2'Ki\Uk-i Uk + i 
fc=i ^ ' fe=i ^ 

In particular, for d ~ I we recover the familiar Szegb kernel S'(z, w) 
^^=zr^ of the upper half plane, and Herglotz kernel becomes 

H (z, u; i) 



ni \u — z u^ + 1 

Let $(2:) be an analytic function, mapping the open upper half-plane into itself 
and satisfying linis^oo ^{is) = 0. Then ln<I>(z) is well defined, with 31n<I>(z) G 
(0,7r). The argument preceding Proposition 7 remains valid, with the result 

^{z) = cxp[iF{z)], 5z > 0, 
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where F{z) is analytic, KF(z) £ (0, tt) and 

F{z) = i'^F{i) + / H{z,u: i)(j){u)du, 
and (j) £ L°°(R), < < tt. In conclusion, we obtain the representation 

a formula already invoked in (21). 

5.2.2. The h ball 

Let /x be a positive Borcl measure supported by the closed unit ball b of M'', 
and denote by 

aaitj) = / x"d/x(x), a g N''. 
We consider a version of the Fantappie transform of /x: 



y. i X ' Z 



where u • v = uivi + ... + UdVd- 

Note that $(/i) is an analytic functions defined in the open unit ball B of C^. 
Its Taylor series expansion at 2; = is reducible, modulo universal constants, to 
the moments of /i: 



Remark also that, for all z G B: 



Jb |1 — X • Z| 

By the maximum principle for pluri-harmonic functions, equality sign can hap- 
pen only if $(/x)(z) is identically equal to a purely imaginary constant, which is 
impossible, since $(^)(0) = /i(b). 
Define the function 

F(z) = -ilni$(^)(z) 

on the unit ball, such that i^{fj.){z) = expiF{z), with ^F{z) € [0,7r] for z a B. 
The function F{z) is analytic in the ball, and has a positive, bounded real part 
there. By the generalized Riesz-Herglotz formula, see [17, 45], we infer: 

F(z) = iSF(0)+ / [2S{z,w)~l](j){w)dcj{w), 

where 9B is the unit sphere, a{w) is the surface element on 9B, normalized to 
have mass equal to one, 

S{z,w) 



(1 - z-wy 
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is Szego's kernel of the ball, and <j>{w) is a measurable function on the sphere, 
satisfying 

< (j){w) < TT, w & dB. 

Let us deal first with the free term, similarly as we did in the one-dimensional 
case: 

$(/.)(0)==Mb), 

and therefore 

F(0) = -iliiifi{h) = -ilnfi{h)+-. 
■^ , ' 2 

iQF(O) 

The total mass of </)(?«) is now easily computed by setting z = in the Riesz- 
Herglotz formula, to obtain 

n f 

— i In t/i(b) + — ■ = — i lni/x(b) + / (f>{w)d(j{w), 

or 

/ (j){w)da{w) = -. 
Job ■^ 

The Fantappie transform of the measure ^ simplified by substituting the 

exact values for the free term and the total mass of the phase function (p: 

i^{fj.){z) = expiF(z) 

= exp ln^(b)+i / [2S{z,w) — l]4}{w)da{w) 

I JdB 

^i^i){z) = /i(b) exp <^ i / [2S{z, w) - l](j){w)da{w) - i- 
[ JdB 2 



= /i(b) exp <.2i [S{z, w) — l](j){w)da{w) 

I JdB 

This formula relates, via a triangular, non-linear transformation, the moments 
(a^(a)) of fi to the moments {a^{a)) of the density (f>. Additionally, note that 
the Koranyi-Pukanszky theorem asserts that all the multivariate moments of 
(j) at mixed-sign indices are zero, i.e., for a niulti- index a, 0^(0;) = 0, except 
possibly when Vi, a^ > 0, or Vi, a^ < 0. 

5.2.3. The h ball 

Let A = {x e R''; |xi| -I- ... -I- \xd\ < 1} and consider a positive measure ^ 
supported by A. Its Fantappie transform is analytic in the unit polydisk D'*, 
and has there a power series expansion 

'a\- 



*(m)W-EM'«a.(")^"- 



Arguing as in the case of the ball, the function ln$(//) is analytic in the open 
polydisk, and has non- negative imaginary part there bounded above by tt. The 
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analog of Riesz-Herglotz formula (as derived for the first time by Koranyi and 
Pukanszky [26]) yields a measurable function (j) defined on the unit torus T'', 
with values in the interval [0, tt], so that: 

$(Ai)(z) =M(A)exp J2i / [W{z,w)~l](t){w)de{w) 



where this time 



n(z,u;) = nY31 



j = l - -'J""! 



and 

den - n '" 



2'iTiwi 
j=l J 



Let a^{a) = /^^ (j){w)w°'d9{w) denote the Fourier coefficients of the function 
(i.e. its trigonometric moments on the torus). 
At the level of generating series we obtain the following transform 



Y, ^°m(«)^" = m(A) exp[2z J2 a0(a)^"]- 



Consider a simple example. 

Example. To verify that the above formulas arc correct, we consider the ID 
case, with the Dirac measure d/i = cSa, where c > 0, a G [—1, 1]. The Fantappie 
transform is the analytic function in the unit disk: 



1 — az 



Since f{z) has positive real part, if{z) has positive imaginary part in the disk. 
Whence lnif{z) is well defined, analytic, and has imaginary part in the interval 
[0,7r], thus the classical Riesz-Herglotz formula is applicable to the function 
~i\nif{z): 

-ilmfiz) = i^[-i\mfm + f i±i^<^(«;)-^. 

Above, 0(u)) — 3ff(— i h\if{w)) is a measurable function on the torus with values 
in [0,7r]. Note that /(O) = c, so that lni/(0) = lnc + i7r/2, and i3[-ilni/(0)] = 
— ilnc. By evaluating z = in the formula we obtain 

— i Inc + 7r/2 = — ilnc + / 4>{'w) , 

Jt 2111111 

that is 

7r/2= / Mw)- . 
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Finally we get 



1 ■ , I-/ \ 1 . f 1 + zw ^ . . dw 

mi + [nj(z) = mc + t / ^<P{'w)- 

Jt 1 — zw 2niw 



ln/(z) = \nc + i I 

Jt 



1 + zw 
1 — zw 



dw 



lnc+22 / ^zH>\w 



In c + 2t 



T 1 — zw ^-KIW 

1 

1 



1 — ZW 



dw 
2TTiw ' 



which is consistent with our general formulas. 

One step further, we can easily compute the positive (n > 0) Fourier coeffi- 
cients oi (t){w) = ^{—i\nif{w)): 

(h(w)w"'- = / --\nif(w)w"-- = 

J ^^ ' 2mw 7t 2 ^ ' 2mw 

If, c dw 1 a" 

In^ w 



2i Jf 1 — aw 2'Kiw 2i n 
The final verification: 



;cxp^ 



oo „ 
n 

n 
Z 



1 — az ^-^ n 

n—l 

5.3. Partial Fantappie transform 

In this section we continue analytically the Fantappie transform of a measure 
supported by a convex cone in a single direction, treating the rest of the variables 
as parameters, with a double benefit: a simple and well known formulas in ID, 
and a tight control of the growth of the phase function. We closely follow below 
the article [23], although similar computations have appeared much earlier in 
the work of Ncvanlinna and Verblunsky. 

5.3.1. Regularization by parametrized single-variable transforms 

The setting is the same: F C M'' is a closed, solid, acute convex cone and 
/x is a finite positive measure supported by its polar cone F*. We consider the 
analytic extension of the Fantappie transform: 



Since 
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Thus, for any fixed y G F, the function z H> i$(— z,y) preserves the upper 
half-plane and hence it can be represented for 3w > as 

where Ci^y) > and < (j)y{t) < 1 both depending measurably on y, respectively 
y and t, see [23]. Both functions C{y),(f)y{t) are uniquely determined by $(z, y), 
hence by the measure fi. For illustration, we provide a simple example of a point 
mass: 

Example. Take n = ccJq, where c > 0. First we obtain by direct integration 

-1 f°° I t 

= exp / (- T— 7j)^^' 

hence 
obtaining 

C(y) = c, 0y(t) = x[o,oo)(0. 2/ e F. 

It is rather annoying that for such a simple measure as a point mass, the phase 
function has an unbounded support. Fortunately, there is a simple remedy, 
derived from an observation of Verblunsky [46, 47]. 

Theorem 9. Let fi be a finite positive measure supported on the cone F*. For 
every y G F there exists a phase function ^y G i^([0, oo),(ii), < ^y < 1, 
measurably depending on y, such that 

^^fMEL=,^^rm±^ 3z>o. (22) 

Jr x-y- z Jq t- z 

Moreover, if j-^* 1^1" dn(x) < oo for some n G N, then /^ t"^y{t)dt < cx) for all 
yGF. 

Proof. Fix a point y G F and denote by fiy the push forward of the measure ^ 
via the map x ^^ x ■ y. Specifically, for a test function / G Co([0, oo)): 



f{x■y)d^i{x)= f{t)d^iy{t). (23) 

r* Jo 

In these terms, Fantappie's transform of the measure /i becomes 

d^iy{t) 



^{z,y) 



z + i 



According to Verblunsky's theorem [23], there exists a function ^y G -^^^([0, oo), di), < 
£,y <1 with the property 



l + $(-z,y)=exp / 
Jo 



^^W^^ Sz>0. 



t- z 
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The boundary limit operations producing ^y from fiy imply that the dependence 
y 1-^ ^y is (weakly) measurable, as a map from F to L-^ {[0 , oo) , dt) . Finally, 
Theorem A.b) of [23] implies the finiteness of the first n + 1 moments of every 
fy, provided that the moments of fiy of the same order are finite. D 

Returning to our simple example, ^ = c6o, we find this time 

^, , _, f cSo(x) ^ c r dt 

1 + $(-z, y) = 1 + / -^^-^ = 1 = cxp ' 



Y* X ■ y — z z Jq t — z 



whence 



Cy == X[o,c], y e r. 

Assume that Jp, |x|" dfi{x) < oo for a positive value of n and denote the 
initial moments of /x as: 

7q = / x°'dii{x), \a\ < n. 



Similarly, denote 



p. 



c(2/)j= / t^^y{t)dt, 0<j<n. 
Jo 



Then one can identify asymptotically the series expansion (in a wedge with 
vertex at 2; = 0) of the two terms in (22), obtaining the algebraic relation 

1-E^(E S^") + 0(.— ) = exp(-X:||?). (24) 

fe=0 |a|=/c i=0 

For details see again [23]. In particular, by equating the coefiiicients of z~^ one 
finds Verblunsky's identity 

co(2/) = 7o, y G r. 

Note that, after expanding the exponential series, the moment Cj{y) is given 
by a universal polynomial function in the variables 7q,, |q;| < j. It is exactly this 
system of polynomial equation which was discovered and exploited by Markov 
(in dimension one). 

5.3.2. Inversion through Radon transform 

Using entropy optimization, for each p £ F, we obtain functions ^* that 
approximate Aronszajn-Donoghue phase functions ^p of push forward measures 
Up defined by (23). The inversion procedure given in Section 4 would be able to 
recover /x* measures, however, piecing together approximation fi* from "slices" 
/x* would be a challenging task. Instead, we seek to recover the Radon transform 
TZfj,*{p,t) of approximation fi, which is then inverted by one of the standard 
algorithms for inverse Radon transform. 
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We start by relating the partial Fantappie transform to Cauchy transform 

$pt(p, ~z) = / ■ = Cfip{z) = expC^p(z) - 1, 

with a slight abuse of notation for the Cauchy transform, see (22). To avoid 
the singular integral in passing from complex to the real Fantappie directly, we 
sum the Plcmelj-Sokhotski formulas (11) to obtain 

"^Kp^Po) = C^p(-po) = hm - [C^ip{-po - ie) + Cfip{-po + ie)] 

ci.0 Z 

= exp [-TTH^pi-po]] cos [TT^pi-po)] - 1, 

using the derivation analogous to derivation of (13), where we defined function 
fp just to relieve notation for the rest of the procedure. 

To connect the Fantappie transform to the Radon transform of a measure, 
we again follow Henkin and Shananin [18]. For a positive measure ^ with density 
<^ supported in the positive orthant M" , define the Radon transform as 

TZfi{0,s) = / ^{x)da{x), (25) 

JHie,s) 

S{s-~x-e)^{x)daix), (26) 



where H{9, s) ~ {x £ R" : x ■ 9 — s} is the integration hyperplane and da{x) its 
surface area measure. Henkin and Shananin give the following relation between 
Fantappie and Radon transforms of rapidly decaying measure fi: 

$M.,Po)=r^^cfr, (27) 

Jo T + po 

valid for po > 0, p € R" . This expression can be interpreted as a composition 
of Hilbert and Radon transforms 

^fi{p,Po)^[nnfi{p,-)]{-Po): 

where Hilbert transform is taken along the offset parameter. As —Ti^ is the 
identity operator, by applying another Hilbert transform, we can evaluate the 
Radon transform as 

'R-KPiPo) = [Hfplipo), 

IT 

where we define 

fpipo) - exp [-1^'Hipijpo)] cos [7r^p(po)] - 1- 

As mentioned before, Hilbert transform can be efficiently evaluated using FFT. 
Therefore, for each selected p, we can evaluate the Radon transform along pq 
axis. One might validly ask why we decided to go through this labyrinthine 
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process only to evaluate a Radon transform of fj,. The answer lies in the method 
used to obtain moment data. If the original measure n has a density fully 
accessible for arbitrary Radon- type measurements, as it is in medical tomogra- 
phy, then the entire procedure is superfluous as we can access 7?.^(p) directly 
at any p € F. However, for singular measures and in some settings, e.g., in- 
variant measures on attractors of dynamical systems. Radon-measurements are 
not directly possible and the described procedure becomes an acceptable path 
to reconstruction. 

An unfortunate obstacle prevents us for completing the process by invoking 
a readily-available inversion algorithm for the Radon transform. The relation 
between the Fantappie and Radon transform (27) holds only in the positive 
orthant {po,p) G R" , which is, in general, not enough for a numerically- 
stable reconstruction [24]. It is possible that this obstacle can be removed by 
considering certain symmetries of the problem. However, these considerations 
would lead us too far from the central theme of this paper and we plan to explore 
them in a subsequent paper. 

6. Conclusions 

In this paper, we presented an approach aimed at representing singular mea- 
sures using bounded densities. Our goal was to use existing entropy optimization 
methods to solve truncated moment problems for singular measures. Previously, 
the entropy optimization could not be applied to singular measures due to lack 
of convergence in the optimization procedure 

The paper adds two steps that bookend the entropy optimization. The regu- 
larization step uses a triangular, recursive transformation on the input moment 
set to produce the conditioned moment set which is a feasible input for entropy 
optimization. The inversion step uses the density that solves the entropy opti- 
mization constrained by conditioned moments, and recovers an approximation 
to the original measure. 

We presented both steps in detail for the support in a one-dimensional space. 
Two different settings, unbounded and bounded supports, were analyzed, re- 
sulting in, respectively, a more general formulation using power moments, and 
a numerically favorable formulation using trigonometric moments. 

We generalize the regularization step to multivariate domains, in particular 
supports of measures in wedges and particular compact domains in R^. The 
inversion, however, proved to be a technically more demanding task, with the 
direct generalization of the Plemelj-Sokhotski formulas requiring a detailed ap- 
plication of theory of Clifford algebras. We did not tackle such generalization 
in this paper, nevertheless, we presented an outline of a simpler tomographic 
process. It would reduce a multivariate inversion problem to a family of inver- 
sion along rays in the domain, which can be solved using the one-dimensional 
method. The information based on ray transforms could then be integrated into 
the approximation of the original measure using methods of tomography. 

This paper is the first part of this research effort. We plan to explore the 
full solution to the multivariate inversion step using tomography in one of the 



33 



follow-up papers. Furthermore, a future paper should formulate error analysis 
and present numerical confirmation of the usefulness of our method, applied to 
concrete examples. 

Appendix A. Miller-Nakos algorithm for exponentiation of a power 



Computing coefficients of an exponentiated (formal) power series can be 
performed recursively, i.e., from knowledge of the coefficients of the base and 
coefficients of the lower coefficients. The original algorithm in a single variable 
is given by Peter Henrici who attributes it to J.C.P. Miller [48]. The extension 
to the multivariate case is due to George Nakos [28], which is unfortunately 
published in a journal that is not easily accessible. Therefore, we give the 
algorithm in its entirety here, stressing that all the original ideas were present 
in the above papers. Ours arc the choice of notation and the (trivial) extension 
to the case of series with zero free terms, mentioned at the end of the section. 

Let a, /3, //, 7 G Ng denote multi-indices and let the basis multi-index be Ci = 
(0, 0, . . . , 0, 1, 0, . . . , 0) where 1 is at the i-th position. Unless noted otherwise, 
sums over multi- indices range over all multi- indices Nq. Furthermore, denote 
d - ^ 

' ~ dxi ■ 

Theorem 10 (J.C.P. Miller, G. Nakos). Let A{x) and B{x) be multivariate 
power series with non-zero free terms, given by expressions 



Aix) 



-E 






B{x) 



= E 

/9 



bjsx^ 



with ao 7^ 0, /3o 7^ 0. 

// the power series are related by equation 

B{x) = A{xf, 

then coefficients bp can be computed recursively by expressions 

bo = a(; 



E- 



0<7<Ai 



'-"^-' 



^7*^/^,-7? 



(A.l) 



where 



x/l3\= Y, "'M 



Proof. The proof of the algorithm rests on the identity 

A{x)diB{x) = kB{x)d,A{x), 



(A.2) 



^Particularly, \a/a\ counts the nonzero elements in a. 
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derived by expanding diB{x) by chain rule and multiplying both sides by A{x). 
Series expansions for derivatives diA(x), and analogously for diB{x), is 

diA{x) — y^aaaiX°'~'^\ 

a 

Expanding (A. 2) into series we obtain 



E«"-") [EhP^^'-A =kr£b,xA (y: 



Ctfj ct^'X' 



Introduce following change of indices: on the left hand side over /3, /3— e^ i— t- /3, 
(3 ^^ j3 + a, Pi 1-^ Pi + 1, and the analogous substitutions on the right hand side 
in sum over a, to obtain the identity 

Y, ac.x'^ I](/3. + l)bp+e,x^' = k^h^^' I](a^ + l)aa+.,x". 

a p pa 

At this point, the goal is to compute coefEcient 6^ from knowledge of Oq and 
h^ for 7 < /I. The products of power scries ^^ a^x'^ ^o hpx^ ~ ^ c^x'^ are 

computed using convolution over coefficients c^ = X]q+/3=7 ^a^/3- The coeffi- 
cients at index fi — Si on both sides of the above identity have to match, yielding 
coefficient identity 

y^ aabp+^.{(3.i + 1) ^ k Y bi3aa+ei{ai + l). 

The expressions can be simplified by another re-indexing, taking into account 
constraints on indices in summation, and symmetry in summation: on LHS, 

a H- 7, Ei + /? ^ ^ — 7, on RHS /3 i-> /i — 7, e^ + a M- 7, to obtain 

or 

y^ a^b^-^ [/ij - (fc + 1)7^] = 0. 

0<7<Ai 

To solve for 5^, extract 7 = case from the sum to obtain 
aab^,^lt = ^^ a^b^^^ [{k + 1)7, - fii] 

0<7</^ 



7 , 0,-fbfj.-^^ 



an 

" 0<7<At 



(fc + 1)^-1 



Since this expression is valid if and only if /i^ 7^ 0, we can sum over all such 
cases to obtain 



Im/mI ^m = — y] a76M-7 [(fc + 1) 17/mI - Im/mI] , 

Ctn — 

0<7<^i 

which, through dividing by |/i//i|, yields the expression (A.l). □ 
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The moment expansion series that we use have a zero free term, which is 
essential for solving for 6^ in the last step of the proof. This can be resolved by 
adding an extra step which uses the binomial expansion 

[(A(.T) + i)-ir= ^ ('"V_i)"-^-[A(x) + if, 

0<fc<n ^'^^ 

where A{x) is the series without a free term. In general, this step does involve 
n extra computations, however, the moment conversions such as (9) require all 
the powers between and n anyway, so there is no additional cost involved for 
our purposes. 



Appendix B. Fast Iterative Algorithm for Entropy Optimization 

This algorithm has been described in [4], and it is based on explicit dis- 
cretization of the entropy functional, which reduces computation of moments 
to a matrix multiplication. For completeness, we present it here in a distilled 
form. 

First, assume the domain is [0, 1] interval, and choose points Xk with quadra- 
ture weights Wk ior k = 1, . . . ,K, i.e., 

1 J 

f{x)dx ^^f{xj)wj. 

For an arbitrary distribution p{x), we will write p = (pj), Pj := pi^j), and 
p = (Pj)i Pj = WjPj. To evaluate moments of p{x) we employ the matrix 
A = {ttij), whose rows are evaluations of monomials on the array Xk- For first 
A'' power moments, A will be a row-truncated Vandermonde matrix a^j = Xp 
where i = 1, . . . , N , and j = 1, . . . , K. The vector of moments n = {pn) for a 
discretized distribution pj is easily evaluated by taking the product // = A • p. 
Let ai,i = 1, . . . , A^ be the set of the Lagrange multipliers (dual variables), in 
which the entropy optimization is unconstrained. The primal is then evaluated 
by function 

pia), = exp [(A^ • a), - l] , (B.l) 

with the goal of finding a* such that pj{a*) « p*. Constraint deviation vector 
is 

h,{a) ^ [A ■ p{a)l^ - fM. (B.2) 

The optimization program strives to achieve 

min (l p(a) — /i a) , 

by cyclically updating components of a. Denote kth iteration of a by a'*"'^. 
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Let fc be a step counter, and set i = [k mod iV) + 1 to be the index of 
the Lagrange multipher updated in the fcth step. Fix the convergence tolerance 
e > 0, and denote initial step by /c = 0. The initial vector a'^^) can be chosen as 
random numbers in some interval, a constant vector, or some other vector. 

Compute the correction factor 

A'^'-ln ., 7,,... , (B.3) 

[A-p(a('^))J. 

and update ith Lagrange multiplier 

a, = a, , for j ^ i. 

If ||/i(a''^+^^)|| < e, then a* = a'^'^+^\ otherwise, increase k by one, and restart 
from computation of the correction factor. 

The paper [4] asserts that the algorithm converges when yu^ > and a.y £ 
[0,1], Vi, or when /i^ < and aij e [—1,0], Vi. When this is not the case, 
the authors provide a pre-conditioning step that modifies A and /i to ensure 
convergence. 

The presented algorithm can be extended to cases where generalized mo- 
ments are taken, instead of power moments. In those cases, the moments are 
not necessarily positive, nor is the matrix A, so they have to be rescalcd before 
the correction factor formula (B.3) can be used. 

To generalized moments on another interval, let Ti{x) be linearly indepen- 
dent functions which are used to generate generalized moments. The matrix A 
is then given by elements Oy ~ Ti{xj) for i~l,...^N,j = l,...^K. Choose a 
positive constant S > and let 

(needed offset) Ui = — min(aij) -t- S 

j 

(scale) Mi ~ max(wj + aij) 

j 

1 



(scaling factor) ti = 



[AU + S) 



The original paper used 6^1 and ti = [N{Mi + S)]^^ but we found that such 
settings result in somewhat slower convergence rates. 

The conditioned matrix A' and moment vector /i' are then computed as 

a,y- = ti (ui + aij ) 
Mj = Uiui + ^J,i), 

which ensures convergence conditions. 

Now the original program is modified by replacing the primal formula (B.l), 
the correction factor formula (B.3), and constraint equation (B.2) by, respec- 
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tively, 

p'(a) = cxp [(A'^ • a), - l] , 

A'(^) = In ^ 

[A'.p'(aW)]/ 

and 

Note that in the constraint deviation, we evaluate the unconditioned moments 
of the conditioned primaL The justification is the fact that p* = p'(a) simul- 
taneously solves the conditioned problem A'p = /i' and unconditioned problem 
Ap = /i. due to linearity of the pre-conditioning step. We, therefore, find the 
unconditioned constraint deviation condition more intuitive to use, which makes 
it easier to choose a desired convergence tolerance e. 
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